Inhomogeneous shear flows in soft jammed materials with tunable attractive forces 
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We perform molecular dynamics simulations to characterize the occurrence of inhomogeneous 
shear flows in soft jammed materials. We use rough walls to impose a simple shear flow and study 
the athermal motion of jammed assemblies of soft particles, both for purely repulsive interactions 
and in the presence of an additional short-range attraction of varying strength. In steady state, 
pronounced flow inhomogeneities emerge for all systems when the shear rate becomes small. De- 
viations from linear flow are stronger in magnitude and become very long-lived when the strength 
of the attraction increases, but differ from permanent shear-bands. Flow inhomogeneities occur 
in a stress window bounded by the dynamic and static yield stress values. Attractive forces en- 
hance the flow heterogeneities because they accelerate stress relaxation, thus effectively moving the 
system closer to the yield stress regime where inhomogeneities are most pronounced. The present 
scenario for understanding the effect of particle adhesion on shear localization, which is based on de- 
tailed molecular dynamics simulations with realistic particle interactions, differs qualitatively from 
previous qualitative explanations and ad-hoc theoretical modelling. 

PACS numbers: 62.20.-x, 83.60.La, 83.80.Iz 



I. INTRODUCTION 

Soft jammed materials, such as dense emulsions, foams 
and pastes are ubiquitous in nature and have a wide range 
of industrial applications 0, Q ■ Normally, these materi- 
als only flow when an externally applied stress exceeds a 
critical value, the "yield stress" , while they behave as a 
soft solid otherwise. Thus, the flow properties of these 
systems are intrinsically non-linear and exhibit complex 
features that challenge both experimentalists and theo- 
reticians [H . To understand the rheology of these com- 
plex fluids, we want to know the mechanical response 
of the system to an externally applied force. However, 
this becomes a highly non-trivial task when the response 
is not spatially homogeneous, or when it occurs over a 
broad range of timescales. These two challenging fac- 
tors routinely characterize the rheology of soft jammed 
materials, emphasizing the need for spatially and tem- 
porally resolved studies of the flow properties in dense 
particulate systems under shear Q. Clearly, molecular 
dynamics simulations are well-suited to pursue this task, 
because they naturally combine particle resolution with 
the possibility to simulate large systems with controlled 
interactions over relatively long times. 

In typical experimental conditions, a broad variety 
of complex fluids display spatially inhomogeneous flows, 
usually described as "shear-bands" , even though this sin- 
gle name in fact hides a diversity of distinct phenom- 
ena [H. In this paper, we are specifically interested in 
the flow behaviour of dense assemblies of large spherical 
particles that form athermal disordered solids, and we 
primarily think of foams, emulsions and dense colloidal 
suspensions as relevant experimental realizations of our 
numerical model Q ■ While early experiments commonly 
reported the existence of inhomogeneous shear-banded 
flows in foams and dense emulsions, more recent work 0- 



has established that, when properly prepared and stud- 
ied over sufficiently long times, shear-bands in dense sys- 
tems of soft repulsive particles do not appear as a per- 
manent phenomenon, although flow inhomogeneities may 
appear to be extremely long-lived in some cases [f| . Re- 
cently, however, it was reported in several instances that 
the addition of a small amount of attractive forces be- 
tween particles triggers the appearance of shear-bands 0- 
[Toj . A common interpretation is that shear- bands in this 
case result from thixotropic behaviour competing with 
the imposed flow [ll( . While natural for low-density col- 
loidal gels which have a complex structure [HI, EH , this 
explanation appears less convincing for jammed systems, 
which present instead a fairly homogeneous structure. 

Shear-banding phenomena have been observed also in 
numerical simulations of amorphous solids under flow (Til . 
Il5j . It was initially reported for a sheared Lennard- Jones 
mixture where it was argued that the coexistence 
of flowing and static phases results from the existence 
of distinct bounds for static and dynamic yield stresses 
leading to a multivalued flow curve |17ll . Subsequently, in 
studies of model metallic glasses [18| , shear-banding was 
observed using a variety of boundary conditions, quench 
rates, or systems sizes, which motivated theoretical ex- 
tensions of the shear-transformation zone model [lj| to 
account for shear-bands (2p| . However, in experiments 
performed on actual metallic glasses, these flow inho- 
mogeneities evolve rapidly with the applied deformation 
and the system develops fractures before a steady state 
can be reached. In a parallel effort, studies of athermal 
quasi-static shear flow of amorphous solids have revealed 
the existence of system-spanning avalanch es g enerated by 
correlated activation of plastic events [U [22( . These ob- 
servations seem to be in tune with the generic scenario 
that dynamical heterogeneities are a characteristic fea- 
ture of amorphous materials 15 1 . Although it is tempting 



2 



to speculate that dynamic heterogeneities, avalanches, 
and shear-bands are various facets of the same underly- 
ing physics, more precise links between these phenomena 
are missing flij ]. 

At the theoretical level, many early models developed 
to account for the rheology of soft amorphous materials 
were mean-field in nature, and spatial fluctuations were 
usually discarded [23T - I25] . More recently, several coarse- 
grained models have emerged that attempt to capture the 
idea, revealed by the above mentioned numerical studies, 
that plastic events are localized but may trigger addi- 
tional plastic events elsewhere in the system, thus cascad- 
ing into sustained flow pril - l30j . While such modelling di- 
rectly yields spatially inhomogeneous dynamics, the ap- 
pearance of permanent shear-bands does not necessarily 
follow. Numerical simulations of these models indeed do 
not produce genuine shear-bands [2(| [I?} , which seems to 
suggest that shear-bands might only occur under quite 
specific conditions. Several recent studies of simple mod- 
els suggest that some form of long-lived shear-banding 
phenomena may occur after shear start-up [2(], [33| . 
These approaches also build on the possibility to observe 
distinct static and dynamic yield stress values, the former 
being enhanced by prolonged aging in thermal glasses. 

Several mechanisms have been put forward to account 
for permanent flow inhomogeneities, which typically re- 
volve around the idea that the stress-strain rate flow 
curve, a = 17(7), is multi- valued. We already mentioned 
the possibility, first discussed in [l7| , that the flow 
curve at finite shear rates is monotonic [for instance of the 
Hcrschcl-Bulkley type with a finite dynamic yield stress, 
<Jd = liiri-y-j-o 17(7)], but that there exists a static branch 
at 7 = extending up to a static yield stress value, <r s , 
larger than the dynamic one, er s > <7d- This opens a 
finite range of stress where the shear rate can take two 
values. Genuine non-monotonic flow curves have recently 
been obtained in various models by including the generic 
idea that yielding dynamics should be self-consistently 
connected to the evolution of the local structure. This 
was done for instance using self-consistent dynamics of 
the effective temperature |31j in the Soft Glassy Rhe- 
ology model [23I, HH, or by incorporating thixotropic 
effects using an additional timescale for structural "re- 
structuration" in schematic [35j or coarse-grained mod- 
els [H, [H, H3] , in which case non-monotonic flow curves 
arise when structural recovery is slower than the charac- 
teristic relaxation time of the system. This hypothesis 
was however not justified by microscopic arguments or 
detailed measurements. Alternatively, it has been pro- 
posed that shear-banding could also occur due to a shear- 
concentration coupling [381 ] . with the possibility that at 
small enough shear rates, variations in local concentra- 
tions would result in the large fluctuations in flow-rates. 
Clearly, more numerical and experimental studies are 
now necessary in order to test and discriminate these 
different physical ideas. 

In this paper, we report simulational studies of the 
athcrmal flow of highly jammed systems, consisting of 



particles having either purely repulsive interactions (as in 
foams [j| and simple emulsions pHH), or having both re- 
pulsive and short-range attractions (as in adhesive emul- 
sions j7l-ficil]). Motivated by the phenomonological find- 
ing that attractive forces might be responsible for the ap- 
pearance of permanent shear-bands, we specifically check 
whether changing interactions results in different shear 
localization properties during steady state flow, and use 
our spatially and temporally resolved simulations to seek 
a physical interpretation of our findings. We find that 
strong flow inhomogeneities are present in all systems, 
except that the lifetime and degree of fluctuations are in- 
deed much higher with increased attractions, but we do 
not observe simple, permanent shear-bands even in our 
most adhesive systems. We find instead that attractive 
forces, by accelerating stress relaxation, effectively shift 
the system closer towards the yield stress regime where 
flow inhomogeneities can become so pronounced that the 
system is not able to sustain a linear flow profile. Such 
a theoretical scenario was was not anticipated in any of 
the above-mentioned work. 

The paper is organized as follows. In Sec. [XT] we 
describe our model and numerical methods. We then 
present our measurements and results in Sec. IIII1 and we 
discuss our results in Sec. [IVJ We conclude the paper in 
Sec. El 



II. MODELS AND NUMERICAL METHODS 

A. Repulsive interactions 

The model system that we study is a collection of 
polydisperse soft particles, introduced as a model for 
foams [3^, H(j , and which has now been extensively stud- 
ied to understand the physics of jam med soft materials [l[ 
both in athermal conditions [4lU42| and at finite temper- 
atures [43| . 

In the repulsive case, two particles, having diameters 
di and dj , interact via a harmonic potential: 




ij ^ dij 



(i) 



where dij ~ (di + dj)/2. We choose a 50:50 binary 
mixture for the polydispcrsity with a mean diameter of 
(d) = 1.0 and a size ratio of 1.4 to avoid crystallization. 



B. Attractive interactions 

We introduce adhesive interactions between the parti- 
cles in a manner similar to models of cohesive granular 
media j44[. Specifically, we introduce two parameters, 
l\ and £2 > t\, through which we can control the range 
and depth of the attractive forces between the particles. 
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FIG. 1: Interparticle potential V(r), Eqs. (HE]), for e = 1/2, 
and tunable strength of attractive forces from it = (purely 
repulsive case) to u = 0.15. The range of the attractive forces 
is constant, r/d — 1.2, but the strength increases with u. In- 
set: The corresponding forces (shown with the same colours) 
with the labelling for u = —l\/d quantifying the strength of 
the attractive part. 



We choose the following form for the interparticle inter- 
actions: 



'I.J 



<l + h 



0, 



(2) 

This simple form is chosen because it yields an interaction 
force which is piecewise linear, see Fig. [1] 

In our simulations, we use £2 to fix the range of the 
interparticle force: £2 = 0.2, and change the depth and 
range of the attractive part by varying l\. The shape of 
the interactions, with the inclusion of the attractive part, 
is illustrated in Fig. [T] By varying the position at dis- 
tance (l+£i)dij of the minimum in the interparticle force, 
the amplitude of the attractive force is varied, which we 
quantified by introducing the parameter u = —2e£i/dij, 
which sets an energy scale for the attraction strength. 
We have studied different values of the particle adhesion, 
u = 0, 0.05, 0.10, and 0.15, as shown in Fig. [TJ We do not 
explore larger adhesion strength to avoid the occurrence 
of any shear-induced phase separation. 



C. Simulation methods 

We study the shear flow of a two dimensional system 
of soft particles. We use very large dimensions for the 
simulation box, L x = 84.46(d) and L y = 99.39(d), which 
contains N = 10404 particles. Shear is imposed via rough 
walls constructed as follows. A layer of particles having 
thickness 2(d) is frozen both at the top and bottom of 
the simulation box in the y direction, from an unsheared 
relaxed configuration. These rough walls have thus a 



structure similar to the sheared system, and the same two 
walls are used throughout this work. In the simulations 
with attractive forces, we use the same structure for the 
walls, but implement the same attractive forces for the 
wall-fluid interactions. 

To impose a constant shear rate 7, we pull the top 
wall at a velocity fixed by f wa ii = i{L y — 4(d)), with the 
bottom wall being kept fixed. We also carry out some 
shear simulations with an imposed constant stress, a. 
This is done by pulling the top wall by a tangential force 
F = <jL x (45l |46| , with the bottom wall again remaining 
fixed. 

The motion of the particles in the bulk arc governed 
by the conservative forces described above, while athcr- 
mal behaviour is ensured using viscous dissipative forces. 
During the flow, when two particles overlap, they expe- 
rience a dissipative force which depends on their relative 
velocity: — b[(Vi — Vj)-rij]fij , where b = 2 is the damping 
coefficient, and r,j is the unit vector between particles i 
and j. The range of the "overlap" used for dissipation is 
dij for pairs of purely repulsive particles, and corresponds 
to (1 + t<i)dij when attractive forces are included. 

We work at a constant volume fraction of <f> = 
Ntt (d 3 ) /(6V) = 1.0, which is much beyond the jamming 
point (f>j ~ 0.85. This implies that the structure is fairly 
homogeneous and resembles the one of dense amorphous 
solids. Thus, our modelling approach bears no similar- 
ities with colloidal gels, which are found to also exhibit 
shear-bands in experimental work fl2l Il3j . 

The units for energy, length and time are 2e, (d) and 
(d) j y^2e/m, respectively, where m is the mass of the 
particles. The trajectories of the particles are evolved 
by numerically integrating the corresponding Newton's 
equations of motion, using a velocity- Verlet scheme 47 1. 



III. MEASUREMENTS AND RESULTS 
A. Repulsive Particles 

As discussed above, several recent experiments have 
suggested that jammed materials made of soft repulsive 
particles do not exhibit steady state shear-bands. For 
the system of harmonic spheres introduced by Durian 
to study wet foams [39l |. simulations have both reported 
the presence [46| or absence [48| of shear localization. 
The same model also displays strong dynamical het- 
erogeneities under shear, similar to unsheared thermal 
glasses, quantified by measuring the usual dynamical sus- 
ceptibility X4 [HI- Such heterogeneities have also been 
recently measured in the flow of NIPAM particles which 
similarly interact solely via repulsive interactions [o0j |. 
Additionally, in the quasi-static limit the same model dis- 
plays system-spanning avalanches [5l| and strongly non- 
affine particle motion (52[. Therefore, there are marked 
spatiotemporal fluctuations and non-affinc particle mo- 
tion for this system under shear flow, which certainly 
need to be taken into account when discussing velocity 
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FIG. 2: Left: Map of non-affine displacements for repulsive 
particles, measured for an imposed shear-rate of j = 2.5 x 
10~ 5 and during a strain window of A7 = 0.10 taken during 
steady state flow. Right: The corresponding velocity profile, 
measured during the same strain window. The dashed line 
shows the velocity profile expected for a linear flow. A clear 
"band" can be observed near the middle of the system, which 
spans the system horizontally, but it is not permanent. 



profiles and their fluctuations. 

We begin to explore kinetic heterogeneities by look- 
ing at the map of non-affine displacements of the parti- 
cles, defined as deviations from the local single-particle 
displacements expected from assuming a linear velocity 
profile. In the left panel of Fig. [3J we show such a map 
corresponding to a measurement done within a strain 
window of A7 = 0.10 during steady state flow at an 
imposed strain rate of 7 = 2.5 x 10~ 5 , while the right 
panel shows the corresponding velocity profile. One can 
clearly observe the spatial heterogeneities in the dynam- 
ics with regions having different mobilities within this 
period of deformation. Moreover, it is evident that the 
particles which have undergone large displacements clus- 
ter together to form a "band" aligned in the flow di- 
rection, similar to what has been for instance observed 
in amorphous Lennard- Jones solids under quasistatic de- 
formation (22|, at finite shear rates at T = [53|, or 
in supercooled liquids HH- Quite often, the appearance 
of such system-spanning heterogeneities are invoked as 
proof of presence of "shear-bands" . 

The spatial variation of mobilities, observed during 
this strain window A7, should also be reflected in the 
velocity profiles measured during the same interval. In- 
deed it does. The velocity profile measured during the 
same strain window when the above map was generated 
is shown in the right panel of Fig. [2] The region of large 
mobility in the map corresponds to a large deviation from 
a linear velocity profile. 
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FIG. 3: A series of representative velocity profiles for repul- 
sive particles, averaged during a strain window of A7 = 0.10, 
and sampled during steady-state flow at an imposed 7 = 
2.5 x 10 -5 . These profiles reveal strong deviations from linear 
profiles, which strongly fluctuate both in space and time. 



A number of velocity profiles, all measured during in- 
dependent strain windows of A7 = 0.10 during steady 
state flow at the same value of the shear rate are shown 
in Fig. [3] Large fluctuations deviating away from the 
expected linear velocity profile are generically observed. 
However, during flow the location of the more mobile re- 
gion is seen to switch from one place to another, and the 
observed "shear-bands" are in fact not permanent but 
have their own dynamics. Moreover, we do not always 
see a clear "coexistence" between two distinct regions of 
mobility, which is often associated with shear-banding. 
Instead, the velocities might sometimes evolve more grad- 
ually across the channel. 

The observation of velocity profiles suggests that be- 
fore drawing conclusions about the presence of shear lo- 
calization it is necessary to study and quantify more pre- 
cisely the degree of fluctuations of the velocity profiles, 
in order to answer the following two questions. How do 
the observed heterogeneities depend on the strain window 
chosen to average the profile? How do these fluctuations 
depend on the imposed macroscopic shear rate? We feel 
that such quantitative information is mandatory when 
reporting on shear-banding phenomena. 

To analyze these fluctuations, we average the velocities 
in the a;-direction to compute the local strain rates j(y) 
averaged over a given strain window A7. One can then 
construct, for any chosen A7, a histogram of these locally 
observed strain rates, which we denote N(j). Clearly 
this probability distribution function depends on the two 
key parameters whose influence we wish to study, namely 
the strain window, A7, and the macroscopically imposed 
shear rate, 7. 

We first show the evolution of N(j) with the strain 
window for a fixed value of the shear rate, 7 = 5 x 10~ 5 
in the top panel of Fig. @] For small strain intervals, 
A7 = 0.04 — 0.10, N(j) spans across a wide range of lo- 
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FIG. 4: Top: Distribution of local strain rates 7(3/) with 
changing the strain window A7 used for the average, at an 
imposed strain rate of 5 x 10~ J . Bottom: Variation of S, 
Eq. Q, which quantifies the spread of the distribution Nft), 
with the strain window A7, for a range of imposed strain 
rates. 



cal strain rates from nearly immobile regions (7 ~ 10~ 8 ) 
to regions which flow faster (7 ~ 3 x 10~ 4 ) than the im- 
posed shear rate. Moreover, the shape of the distribution 
is clearly not symmetric and non-Gaussian, with the ap- 
pearance of a flat tail towards small values of the shear 
rate. This behaviour is clearly consistent with our obser- 
vation that spatial fluctuations in the profiles, shown in 
Fig. [21 actually span the entire system. However, the dis- 
tribution narrows down and becomes closer to a Gaussian 
with an increase of the strain window over which velocity 
profiles are averaged, suggesting that spatial fluctuations 
become less correlated over time. For the largest defor- 
mation, A7 = 2, the fluctuations around the imposed 
shear rate have become quite small. Our observations 
are similar to what was reported for a sheared Lennard 
Jones glass [55[. For small times, only a few plastic events 
occur resulting in the initial large heterogeneities which 
are localized in space. However, if one waits long enough, 
the plastic events proliferate across the system and the 
heterogeneities therefore are erased, which eventually re- 
sults in homogeneous flow. Intriguingly, such description 



of the transient character of shear inhomogencities is also 
reminiscent of the temporal evolution of kinetic hetero- 
geneities characteri zing the structural relaxation of ther- 
mal glassy systems [15( . 

These results suggest that some form of "shear- 
banding" exists in the present system, but flow localiza- 
tion is not a permanent phenomenon. It is thus natural 
to ask about the lifetime of these inhomogencities. To 
this end, we introduce the "dispersity" 5 of the distribu- 
tion AT (7), as the ratio of the standard deviation to the 
mean of the distribution: 



S = 



(7) 



N 



(3) 



where the average (■ • ■ ) at is performed over the probabil- 
ity distribution of the local shear rate, AT (7). The dis- 
persity S is the most natural way to quantify the width 
of this distribution. 

In the bottom panel of Fig. [H we show the variation 
of S with the strain interval A7 for a range of imposed 
shear rates from 7 = 2.5 x 10~ 5 to 10~ 3 . Following from 
our discussion above, we find that S decreases with A7 
for all 7. Therefore, if one averages the velocities for long 
enough strain windows, say larger than A7 ~ 2, nearly 
linear velocity profiles will be observed in all cases. In 
fact, it is interesting to note that for all values of 7, het- 
erogeneities become negligible at approximately the same 
strain interval A7. However, in the regime of smaller 
A7, we see that flows become more heterogeneous with 
decreasing shear rates. 

From these results, we conclude that our model of a 
jammed system with soft repulsive interactions does not 
produce permanent shear localization, although strong 
flow inhomogeneities are detected when insufficient aver- 
aging of the velocity profiles is performed, which might 
be a relevant issue in experiments. This analysis also sug- 
gests that a discussion of shear-banding in soft jammed 
materials can not be separated from a discussion of their 
spatio-temporal dynamics. In particular, we have pre- 
sented in Fig. [4] a simple method to quantify the lifetime 
of these inhomogeneities. In Sec. IIVI we will relate these 
observations to the global rheology of the system. 



B. Including attractive interactions 

Motivated by recent experimental results @, [1] where 
the inclusion of particle adhesion in dense emulsions re- 
sulted in qualitatively different flow patterns compared 
to repulsive emulsions, we proceeded to explore the na- 
ture of flow heterogeneities in sheared soft disks with the 
tunable attractive interactions shown in Fig. [TJ 

For these sticky particles, we again look at the map of 
non-affine displacements. In the left panel of Fig. [5j we 
show the spatial map of such displacements for u = 0.15 
during a strain interval of A7 = 0.10, measured in steady 
state at an imposed shear-rate of 7 — 10 -4 . The corre- 
sponding velocity profile is shown in the right panel of 



FIG. 5: Left: Map of non-affine displacements for attractive 
particles, u = 0.15, measured for an imposed shear-rate of 
7 = 10~ 4 and during a strain window of A7 = 0.10 taken 
during steady state flow. Right: The corresponding velocity 
profile, measured during the same strain window. The dashed 
line shows the velocity profile expected for a linear flow. Com- 
pared with Fig. [21 the flow is much more inhomogeneous, the 
top 80 % of the system being nearly unsheared, while the 
bottom part near the wall is sheared very strongly. 

Fig. [5] In the bottom of the map, we can again very 
clearly see a "band" formed by the most mobile particles, 
which spans across the entire length of the system in the 
flow direction, and having a transverse width of around 
10 — 20 diameters, while the top of the system is mainly 
unsheared. The corresponding velocity profile also re- 
flects this via its strongly nonlinear shape in the entire 
system. More interestingly, the mobile band populates 
the bottom of the shear cell adjacent to the static wall, 
whereas the quiet ones are adjacent to the top wall of the 
cell via which shear is generated across the cell, which is 
also evident from the corresponding vclocty profile. 

This is further illustrated by looking at a set of con- 
secutive velocity profiles measured during the same strain 
window, A7 = 0.10, as shown in the top panel of Fig. |6] 
We observe quite dramatic deviations from a linear ve- 
locity profile in all cases and one can clearly see a switch- 
ing of the position of the more mobile population with 
time as the shear-band continuously flips from one wall 
to the other. Note that the shear rate at the center of the 
channel is always very small, such that the bulk of the 
system either flows with the right wall, or remains immo- 
bile with the left wall. Although seemingly reminiscent 
of the velocity oscillations reported for colloidal particles 
in microchannels [5(| , the motion we observe with attrac- 
tive particles is instead very far from periodic, and it is 
enhanced by low rather than fast shear rates. 

When averaged over a much longer period, A7 = 4, 
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FIG. 6: A series of representative velocity profiles for attrac- 
tive particles with it = 0.15, averaged during a strain window 
of A7 = 0.10 (top), and A7 = 4.0 (bottom), and sampled 
during steady-state flow at an imposed shear rate 7 = 10~ 4 . 
These profiles reveal strong deviations from linear profiles, 
which strongly fluctuate both in space and time. Compared 
to the repulsive system, these profiles remain strongly non- 
linear at large deformation A7 = 4.0, suggesting that a linear 
velocity profile is not stable in the presence of strong particle 
adhesion. 



the velocity profiles still deviate significantly from linear 
profiles, as shown in the bottom panel of Fig. [5] For the 
profiles measured during this period, the central region 
appears nearly unsheared, whereas the regions close to 
both the top and bottom walls show local shear rates 
which arc larger than the imposed value. Thus, the flow 
heterogeneities arc much more pronounced in the attrac- 
tive system, with much clearer signs of the "coexistence" 
between sheared and unsheared regions, and this inhomo- 
geneity seems to be persistent over much longer strain 
intervals for the attractive systems as compared to the 
repulsive one. 

Are these "shear-bands" permanent objects? To an- 
swer this question for these sticky particles, we again 
characterize their lifetime, repeating the exercise per- 
formed for the repulsive particles. For each strength 
of attraction it, we compute local strain rates 7 from 
velocity profiles averaged over different strain intervals 
A7. We build the corresponding histograms, Nft), and 
measure the dispersity, 5, of the local strain rates from 
Eq. ([3]). Since the lifetime of the inhomogeneities be- 
comes large when attraction increases, it becomes numer- 
ically difficult to sample a large number of independent 
fluctuations and get accurate statistics for the distribu- 
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FIG. 7: Effect of particle attraction on the variation of the 
dispersity 5, Eq. (j3]), which quantifies the spread of the dis- 
tribution -/V(-y), with the strain window A7, for an imposed 
shear rate of 7 = 1CP 4 . Attraction enhances both the ampli- 
tude and the lifetime of the flow heterogeneities. 
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FIG. 8: Global flow curves a = (7(7) for the repulsive parti- 
cles with u — 0.0 (left) and attractive particles with u = 0.15 
(right). The data are shown with symbols, while the dashed 
lines are fits to the Herschel-Bulkley form, Eq. Q. Horizon- 
tal arrows indicate the stress values at which constant stress 
simulations are performed, see Figs.[!|]and[TU] Notice the very 
different stress scales used in both panels. 



IV. DISCUSSION AND INTERPRETATION 



tions. This forces us to impose a larger shear rate, 7 = 
10 -4 , because much longer simulations would be needed 
to obtain good statistics at lower shear rates, where, 
presumably, even longer-lived and more pronounced are 
present. 

The data for 7 = are shown in Fig. [7J We can 

distinctly see that for all strain intervals, increasing at- 
traction between the particles results in increased hetero- 
geneities. We suspect the effect would be even more pro- 
nounced for lower shear rates. In fact, within the strain 
windows that we have been able to sample, the hetero- 
geneities have not died out when u is large enough. For 
instance, the dispersity obtained over an averaging win- 
dow A7 = 2.0 for 11 > 0.12 is larger than the dispersity 
for u = and A7 = 0.1. However, 5 is still clearly a 
decreasing function of A7, which suggests that for even 
larger strain windows, the temporal fluctuations reported 
in Fig.[S]eventually decrease the overall width of the shear 
rate distribution N(j). While wc conclude that there 
is no permament shear-bands in our jammed adhesive 
systems, our data are nevertheless in broad agreement 
with experimental results showing that particle adhesion 
strongly enhance localization of the shear @-Q3, in the 
sense that our most adhesive systems do not sustain sta- 
ble linear flow profiles, even when averaging velocities 
over deformations as large as 400 %. 

In the following section, wc shall discuss the physical 
origin of these observations. 



A. Static and dynamic yield stresses 

To understand the above observations and the effect of 
the attractive forces, we now turn to the global flow curve 
of the system and ask which (if any) of the theoretical 
arguments summarized in the introduction might apply 
to our system. 

In Fig. we show the flow curves a = (7(7) obtained 
by averaging the stress over the entire system and very 
long times in our simulations with various values of the 
imposed shear rates. For both repulsive (u = 0) and 
strongly attractive (u = 0.15) interactions, we find that 
the resulting flow curves are monotonic functions, thus 
ruling out the possibility that non-monotonic flow curves 
could arise in dense athermal systems with or without ad- 
hesive forces. Additionally, we find that all flow curves 
are well-described by an Herschcl-Bulklcy fitting func- 
tion, which wc write as: 

tfhb = o d (1 + (r c7 )") , (4) 

which contains three fitting parameters. The first pa- 
rameter is the dynamic yield stress, <7d, defined as the 
7 — > limit of the measured stress. The second parame- 
ter, n, determines the shear-thinning behaviour observed 
for shear rates that are large enough to be away from the 
"yield regime" (defined by a pa Ud). The third param- 
eter, t c , has the dimensions of a timescale. It indicates 
at which imposed shear rate the flow is close to the yield 
regime, 777 <C 1, or far from it, 777 3> 1. The obtained 
values of these fitting parameters arc further discussed in 
Sec. lTVBl 
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In the study of the formation of shear-bands in 
Lennard- Jones glass [l6|, evidence was found of a multi- 
valued flow curve by including also the behaviour of the 
system at rest, 7 = 0. This multi-valued nature was due 
to the fact that the value of the static yield stress, a s , 
determined by checking for the onset of flow by applying 
an increasing external stress on a quiescent amorphous 
state [HI, was found to be higher than the dynamical 
yield stress ad defined from finite shear rates measure- 
ments. A strict inequality, a s > cr^, indeed opens a stress 
window, a £ [ud,<J s ], where the system can either be at 
rest, 7 = 0, or flow at finite rate 7 > 0, with the possi- 
bility that both solutions coexist in space, thus possibly 
giving rise to flow inhomogeneities [171 ]. 

Having ruled out non-monotonic flow curves at finite 
shear rates, we thus decided to check whether this sce- 
nario is a possible explanation of the observed flow het- 
erogeneities. Unfortunately, determining a s numerically 
is not an easy task (45j . To proceed, we performed simu- 
lations with an imposed constant stress, a, to determine 
whether the static and dynamic yield stresses differ for 
our systems. By definition, the system should not flow 
when a < a s . 

Using the global flow curves shown in Fig. [51 we se- 
lect for it = two different values of stress to carry 
out our constant stress simulations. The first value is 
a = 0.00518, which corresponds to 7 = 2.5 x 10~ 5 . 
Second, we use a = 0.00566, which corresponds to 
7 = 10~ 4 . Both values arc indicated with horizontal 
arrows in Fig. [8j and the results for these runs are pre- 
sented in Fig. [H To determine whether the system flows 
or not, we simply measure the average velocity Wbuik of 
the system in the direction of the applied stress. As ini- 
tial conditions for these runs, we choose arrested states 
corresponding to local minima in the energy landscape 
of the soft disks. 

The top panel in Fig. [9] shows that for all the different 
trajectories at a = 0.00518, the motion of the fluid par- 
ticles quantified by the average velocity ^buik eventually 
stops at long times. On the contrary, for a = 0.00566, 
the flow continues at some finite velocity Wbuik- Thus, this 
tells us that a s lies in between these two stress values, 
and that indeed it is larger than the dynamic yield stress 
(Td for this system. The situation is therefore similar to 
the observations reported for a Lennard- Jones glass [l(| . 
The fact that there is no flow at a = 0.00518 indicates 
why the heterogeneities are more significant during the 
simulations at 7 = 2.5 x 10~ 5 , in comparison to the flow 
at 7 = 10~ 4 , see Fig.H 

However, one should note that earlier numerical stud- 
ies have reported that the difference between a s and ad 
seems to decrease (albeit quite slowly) as the system size 
is increased towards the thermodynamic limit |46l 
Thus, this indicates that the corresponding degree of flow 
heterogeneties could also have finite size effects, despite 
the fact that our systems already comprise a large num- 
ber of particles, N w 10 4 , and that the channel we use 
is about 100 particles wide. We also note the following 
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FIG. 9: For repulsive particles, u = 0, velocity rwk of the 
center of mass of the system in the direction of the imposed 
shear stress starting from a configuration at rest. Different 
colors correspond to independent initial conditions. Compari- 
son between both panels shows that the static yield stress lies 
in between the two simulated stress values, a = 0.00518 < 
a B < a = 0.00566. 



paradox. Given that the lengthscale of kinetic hetero- 
geneities grows as the shear rate decreases and is believed 
to diverge as 7 — > in systems characterized by a yield 
stress [49(, for any finite size system there should exist a 
shear rate below which finite size effects become relevant, 
and the mechanism mentioned above for the appearance 
of strong flow inhomogeneities could then become rele- 
vant as well. 

We now switch to the attractive system with u = 0.15 
and ask whether forcing the system with a stress value 
corresponding to 7 = 10~ 4 induces flow in the system as 
it does for repulsive particles. The corresponding mag- 
nitude of the external stress, a = 0.04858, is indicated 
by an arrow in the right panel of Fig. [8J The resulting 
data are plotted in Fig. [TU1 and should be compared to 
the bottom panel of Fig. [HI We observe that for all ini- 
tial states, the system eventually comes to rest. Thus, 
in this case, the static yield stress a s is larger than the 
shear stress corresponding to 7 = 10~ 4 . This result is 
in agreement with our earlier observation that at this 
shear rate, attractive systems display much stronger flow 
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FIG. 10: Same as Fig. [9] for attractive particles, u — 0.15, 
and a stress value shown with a horizontal arrow in Fig. [5] 
corresponding to 7 = 10~ 4 . For all initial states, the system 
remains at rest, showing that a s > 0.04858. 



inhomogeneities than repulsive ones, recall Fig. [7J 

Therefore, we conclude from this section that the ex- 
istence of distinct values for static and dynamic yield 
stresses accounts well for flow inhomogeneities in our sys- 
tem, as opposed to a non-monotonic flow curve at finite 
shear rates. Moreover, we also showed that the stress 
window where this competition becomes relevant corre- 
sponds to shear rates values that become larger when 
attraction is increased, as shown by our constant stress 
simulations. We finally remark that these latter results 
suggest that by making simulations at constant shear 
stress, the fluid velocity would remain zero as long as 
a < <t s , and would jump discontinuously to a finite value 
at a s . This is nothing but the "viscosity bifurcation" ob- 
served in several experiments d, [H, [59| . In this language, 
our results imply that increasing the attraction increases 
the value of the "critical" shear rate 7 C at which a steady 
state flow appears, which seems consistent with experi- 
mental results 

Thus, we come to the conclusion that the effect of 
particle adhesion is to continuously increase the critical 
shear rate above which stable linear profiles exist. There- 
fore, for adhesive particles, the regime below the critical 
shear rate becomes easily accessible in numerical simula- 
tions, leading to observation of more prominent inhomo- 
geneities in flow. However, we do not obtain numerical 
evidence that some novel physics comes in, such as for 
instance a slower restructuration process, as advocated 
in Refs. JH-H^). 

We now seek a more microscopic explanation of these 
effects. 



B. "Universal" rescaling of flow curves 

We thus see the emergence of a universal scenario 
for the existence of flow heterogeneities, irrespective 
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FIG. 11: Scaled representation of the flow curves for all pa- 
rameters, using the scaled variables in Eq. ([5}. The dashed 
line is the simple functional form in Eq. ©. Inset: Variation 
of the dynamic yield stress ad (left), and of the timescale r c 
(right) with the strength of the attractive forces u. 



of the interparticlc interactions. Firstly, we observed 
that for both attractive or repulsive interactions, the 
heterogeneous flow corresponds to the regime bounded 
by values of static and dynamic yield stresses. More- 
over, the flow curves for all the different systems (u = 
0.0,0.05,0.10,0.15) are well fitted with the same func- 
tional form, Eq. We find empirically that the shear- 
thinning exponent n varies very little around the value 
n w 0.5 with no systematic trend. Thus, we fix n = \ in 
the following and determine <7d and r c from a fit to the 
data. This analysis thus suggests that the flow curves 
for all our systems can be collapsed by using the scaled 
variables 



a 

y = — ; x = 

onto the simple functional form: 



y = 1 + sfx. 



(5) 



(0) 



This procedure is applied in Fig.[TTJ where the two insets 
display the evolution with the strength of the attraction 
of the two parameters and r c obtained by fitting the 
flow curves. 

We find that the yield stress 072 increases with the at- 
traction u, which is expected as the attraction should 
indeed make the emulsion more cohesive and harder to 
deform. A more surprising result is that the timescale r c 
is found to decrease dramatically by about two orders of 
magnitude between u = and u = 0.15. We shall dwell 
on this unexpected result in the following section. 

The significant observation following the data collapse 
shown in Fig. [IT] is that, for the range of imposed shear 
rates explored in our simulations, the data points corre- 
sponding to increasing attraction correspond to dramat- 
ically decreasing values of the scaled variable x — T c j in 
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the master curve. Thus, varying the degree of attraction 
allows us to "slide" along this master curve, so that the 
most attractive system is effectively much closer to the 
yield stress regime, where strongly heterogeneous flow 
can be expected. On the other hand, repulsive particles 
correspond to a more fluidized segment of the flow curve, 
which explains why only relatively short-lived inhomo- 
geneities are observed in that case. 

While framing their model for shear-banding in 
Ref. [35|, the authors obtained a similar rescaling of ex- 
perimental data for varying degrees of loaded emulsions. 
However, they attribute an increasing tendency towards 
shear-banding to an increasing timescale for "restruc- 
turation" of the material under shear with higher degree 
of loading (stickiness), which is a trend opposite to what 
we observe here. For clarity, we should emphasize again 
that at the large value of the volume fraction (<f> = 1.0) 
that we have studied here, the material is not a gel but 
a highly homogeneous material. For this reason, we be- 
lieve that local fluctuations in the volume fraction are not 
significant, which also rules out the flow-concentration 
coupling scenario proposed in another model [38j for our 
jammed systems. 

The above analysis implies that varying the attrac- 
tion changes the nature of inhomogeneities by affecting 
the scaled variable x — T c j through the evolution of r c . 
Clearly, the same result would be obtained by varying 
instead the shear rate for a given value of t c . This rea- 
soning is consistent with the data shown in Fig. U which 
established that inhomogeneities increase when the shear 
rate decreases. 



C. Microscopic interpretation of timescale r c 

The final piece of information that we need to gather is 
a microscopic understanding of the strong variation with 
the strength of attractive forces of the intrinsic timescale 
t c appearing in the global flow curves, Eq. ((4]). Under an 
applied shear, the system constantly renews its structure. 
We believe that r c quantifies the timescale needed for the 
local stress relaxation to occur. Thus, when the typical 
timescale 1/7 associated with the imposed shear flow is 
much slower than t c , i.e. when t c 7 <C 1, the system is 
effectively in the near quasi-static regime, and a k, a d . In 
the opposite regime, T c j ^ 1, the system has not enough 
time to relax the stress whose averaged value increases 
above the yield stress. 

In this view, r c is a timescale associated to relaxation 
of the stress after deformation. To confirm this inter- 
pretation, we conduct the following numerical experi- 
ment. For the different values of attraction from u = 
to u = 0.15, we apply an instantaneous external strain 
in the x-direction of magnitude e, starting from an ini- 
tially relaxed state (which corresponds to a local minima 
in the energy landscape). We then allow this strained 
configuration to relax [60. [6l|. To fix the value of e, we 
seek a compromise between a very small value where only 




FIG. 12: Relaxation of the stress after a step strain of ampli- 
tude e = 0.2 for the different strengths of interaction u. The 
amount of relaxed stress and the timescale for relaxation both 
decrease strongly with increasing the attraction u. 

trivial elastic deformation would occur, and a very large 
value which would amount to starting from a fully ran- 
dom configuration. We have studied two values, e = 0.1 
and e = 0.2, because they typically correspond to the 
magnitude of the strain whe re p lastic deformations occur 
in quasi-static simulations [51] . We find similar results 
in both cases and show results obtained for e = 0.2. 

We present our data for stress relaxation after these 
step strains in Fig. 1121 where each curve is averaged over 
several (typically 10) initial states. To ease the compar- 
ison betwen different values of the attraction, we plot 
a(t)/a(0), where <r(0) is the stress value right after the 
step strain has been imposed. 

Several remarks can be made based on this figure. Al- 
together the stress exhibits a first rapid decay, followed by 
a plateau at long time scale. First, the short timescale for 
the stress relaxation is observed to decrease strongly with 
increasing the attraction strength. This is in line with the 
result for the time-scale r c in Fig llll Second, it is evi- 
dent that the degree of stickiness influences the amount 
of stress that is relaxed after the step strain, which be- 
comes smaller for stronger adhesion: stickiness increases 
the amount of residual stress that the system can store. 
A possible interpretation is that after deformation ad- 
hesive particles very rapidly "stick" to the neighboring 
particles to minimize locally the potential energy, and 
remain subsequently in this local minimum, while more 
collective, slower moves relaxing larger amounts of stress, 
are more likely to occur for purely repulsive particles. 

In the present system, shear localization is enhanced 
by adhesion because attractive forces accelerate the 
timescale for restructuration after deformation, thus ef- 
fectively moving the system closer to the quasi-static, 
yield stress regime where heterogeneities are more pro- 
nounced. Furthermore the increase of residual stresses 
for adhesive particles does stabilize the transient shear- 
bands over longer times. 
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V. CONCLUSION 

We carried out simulations of two-dimensional jammed 
systems to study the nature of shear flow heterogeneities 
and the influence of the degree of attractive interaction 
between particles. We demonstrated that, independent 
of the nature of interactions, the flow behaviour could 
be described by a universal flow curve, with increasing 
attraction resulting in a flow which is more and more in- 
fluenced by the proximity to the yield stress regime. By 
entering a regime of stress values bounded by the static 
and dynamic yield stress, shear localization is strongly 
promoted, and with increasing attraction, one observes 
enhanced and very long-lived flow inhomogencities. For 
the most strongly adhesive particles, we find that veloc- 
ity profiles do not become linear even when averaged over 
very large deformations, suggesting that in these systems 
a linear flow profile is actually unstable. However, these 
inhomogencities do not take the form of simple, perma- 
nent shear-bands. 

Our results are reminiscent of long-lived flow inhomo- 
geneities measured experimentally in simple yield stress 
fluids @. Interestingly, similar to our adhesive soft 
system, carbopol is measured to exhibit a monotonic 
Herschel-Bulkley flow-curve @, and does not exhibit 
thixotropy, nor aging. An important difference is that 
the experiments report long transients when the shear is 
turned on, while we worked here in steady state. On the 
theoretical side, similar conclusions about the transient 
nature of shear-bands were reached in Refs. pol . lo2l H3| 
on the basis of spatially resolved coarse-grained models 
using the idea that after long aging the static yield stress 
value can become appreciably larger than the dynamic 
yield stress. We emphasize that, by contrast to these 
works, the distinct values for these two yield stresses we 
report here do not depend on the aging time since our 
systems are fully athermal. 

By contrast, fully permanent flow localization tak- 
ing the form of simple shear-bands has been predicted 
through several mesoscopic models [HI, HH, lo5l - l37j [ , which 
all attempt to describe a self-consistent coupling of the 



yielding mechanism to the structural reorganisation {e.g., 
softening mechanisms, timescales separation in the struc- 
tural relaxation, etc.). While these descriptions lead to 
the occurence of permanent bands, this flow behavior is 
systematically associated with the appearance of a non- 
monotonic behavior of the flow curve 17(7) at finite shear 
rate. Although experimental results, simulations and 
the various theoretical models all suggest that this non- 
monotonicity is a necessary condition to observe perma- 
nent shear-bands, we demonstrate here that very long- 
lived, strongly non-linear flow profiles can be observed 
without it. 

Clearly, to make progress and to go beyond the 
above observations and conflicting predictions, it would 
be desirable to make explicit connections between 
the microscopic physical parameters in experiments or 
simulations-such as the adhesive forces considered in this 
work-, and the mesoscopic phenomcnological quantities 
introduced in the various simplified theoretical models, 
such as for instance the distinct timescales which enter 
the definition of coarse-grained elasto-plastic models [35l — 
H3, or the effective 'noise' temperature in the SGR mod- 
els and its numerous variants [3l|, HH . 

Altogether, our results suggest that "shear localiza- 
tion" actually denotes a broad variety of physical be- 
haviours, and further experiments performed with con- 
trolled model systems with tunable adhesion would be 
required to investigate the intimate connection between 
structural recovery and shear localization along the lines 
of the present work. Similarly, it would be interesting 
to design a simple microsopic model yielding the type 
of permanent shear-bands envisioned by coarse-grained 
models, for instance by studying simple atomistic mod- 
els for colloidal gels. 
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